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Abstract 

Most numerical solvers and libraries nowadays are implemented to 
use mathematical models created with language-specific built-in data 
types (e.g. real in Fortran or double in C) and their respective ele¬ 
mentary algebra implementations. However, built-in elementary alge¬ 
bra typically has limited functionality and often restricts flexibility of 
mathematical models and analysis types that can be applied to those 
models. To overcome this limitation, a number of domain-specific 
languages with more feature-rich built-in data types have been pro¬ 
posed. In this paper, we argue that if numerical libraries and solvers 
are designed to use abstract elementary algebra rather than language- 
specihc built-in algebra, modern mainstream languages can be as ef¬ 
fective as any domain-specific language. We illustrate our ideas using 
the example of sparse Jacobian matrix computation. We implement 
an automatic differentiation method that takes advantage of sparse 
system structures and is straightforward to parallelize in MPI setting. 
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Furthermore, we show that the computational cost scales linearly with 
the size of the system. 


1 Introduction 

Numerical analysis is the backbone of engineering design and simulation 
tools today. Solving systems of nonlinear equations, ordinary differential 
equations, or differential-algebraic equations numerically is the key part of 
model-based design. In many instances, the model developer needs to sup¬ 
ply system Jacobian and sparsity pattern to the numerical solver in addition 
to model equations. Since all the information required to identify system 
connectivity structure and assemble the Jacobian is contained within model 
equations, a number of attempts has been made to automate Jacobian gen¬ 
eration. 

Perhaps the most commonly used library for automatic differentiation is 
ADOL-C [1], currently developed at University of Paderborn. ADOL-C is 
built with efficiency in mind and works well with legacy numerical solvers. 
The library includes drivers for Jacobian and Hessian computation. However, 
it does not support MPI-based parallel execution. Sacado [2] is an automatic 
differentiation package which is part of the Trilinos library [3]. It uses ab¬ 
stract elementary algebra to implement automatic differentiation, but it still 
does not support sparse derivatives. It has to allocate memory for derivatives 
with respect to all system variables. 

In this paper, we present an approach for sparse automatic differentia¬ 
tion which can be parallelized in a straightforward fashion when using MPI 
framework. Our approach is based on dehning abstract linear algebra, sim¬ 
ilar to the approach used in Sacado. We provide prototype implementation 
in C-|--|- and demonstrate linear scaling of the computational cost with the 
problem size, in the serial and parallel case. Furthermore, we show that our 
approach allows for model reconhguration at runtime and overall better code 
reuse in scientihc applications. 

The paper is organized as follows: In Section 2 we discuss motivation for 
looking into abstract elementary algebra interfaces. In Section 3 we describe 
our method for sparse automatic differentiation and we outline prototype 
implementation in Section 4. Preliminary benchmarking results are presented 
in Section 5. Future research directions are discussed in Section 6. 
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2 Motivation 


There is a common misconception within the engineering community that 
precompiled numerical models are “hard-wired” models. This was certainly 
the case (quite) a few years ago when the models were coded in languages 
with few object-oriented features, such as Fortran 77 or C. Many legacy nu¬ 
merical solvers indeed require that model structure is known at compile time. 
This is in part due to legacy solvers using built-in data types to represent 
model parameters and variables and by that also inheriting all the limita¬ 
tions of those data types. In many engineering applications, especially in 
design problems, model structure can (and often does) change at runtime. 
For example, when designing a heat exchanger one may want to keep inlet 
and outlet temperatures constant at operating conditions and optimize for 
heat exchanger geometry parameters. In transient simulations, however, heat 
exchanger geometry is fixed and temperatures are system variables. Variable 
and parameter designation is selected by the designer as needed at runtime. 
Having the ability to reuse the same model for different types of analyses 
is extremely important for streamlining engineering processes and reducing 
model verihcation and validation efforts. 

To address this problem, domain-specihc modeling languages based on 
symbolic code manipulations [4, 5] have been introduced. Tools built around 
these languages allow engineers to work in a more interactive design environ¬ 
ment where they can make modifications of their models at runtime. Under 
the hood, the model encoded in the domain-specihc language is processed 
symbolically and code compatible with the numerical solver is generated and 
compiled on the hy. Then, such hard-wired precompiled model is simulated 
and the result is returned to the user. Symbolic manipulations and compiling 
automatically generated code on the hy allow one to reuse models coded in a 
domain-specihc language for diherent types of simulations and analyses using 
pretty much any kind of numerical solver. The downside of such approach 
is that one needs to support a whole new language that may be quite com¬ 
plex and that the model needs to be regenerated and recompiled every time 
the model structure is modihed. Since the language is domain-specihc, the 
user base is relatively small and there are fewer shared resources available. 
Scaling up this approach to more complex problems is another challenge as 
symbolic preprocessing of model equations may become a bottleneck. In a 
framework such as this one needs to support two diherent parallelization 
schemes - one for symbolic preprocessing of model equations and another 
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one for numerical solving of those equations. Symbolic transformations are 
generally nontrivial to parallelize. The more features the domain-specihc lan¬ 
guage offers, the more complex symbolic processing algorithms become and 
so does their parallel implementation. At the time of this writing, we are not 
aware of MPI-based parallel schemes for symbolic processing of mathematical 
equations. 

Modern object-oriented languages such as C-I--I- and Java, which sup¬ 
port operator overloading, template specialization, type traits, and other ad¬ 
vanced features allow one to create numerical models that can be reconhgured 
at runtime. The same functionality provided by symbolic preprocessing of 
model equations can be implemented by creating custom data types and ap¬ 
propriate libraries in mainstream object-oriented languages. Recently, solver 
frameworks that use abstract data types were proposed [6]. Those frame¬ 
works do not require specihc data types to be used, but only specify elemen¬ 
tary algebra that the data types have to support. By designing models and 
solvers to use abstract data types, one can reuse the same models and solvers 
for multiple analysis types such as forward simulations, optimization, sensi¬ 
tivity analysis, or embedded uncertainty quantihcation. Switching between 
these may be accomplished simply by changing the data type (or configu¬ 
ration of the data type). Furthermore, abstract data types can be used to 
compute automatically the system connectivity graph which could then be 
utilized, for instance, to partition the system into smaller subsystems, per¬ 
form index reduction for differential-algebraic equations, implement tearing 
algorithms, and many other calculations. 

In the remainder of the paper, we focus on automatic differentiation im¬ 
plemented using abstract data types with the application towards automatic 
generation of sparse Jacobians for steady state and dynamic simulations, as 
well as optimization. 

2.1 Problem Description 

Many commonly used numerical libraries for solving systems of nonlinear 
equations, differential equations, or optimization problems require the user 
to provide a Jacobian matrix in addition to model governing equations. Given 
a nonlinear problem of the form 


f{x;p) = 0, 


( 1 ) 
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where a: is a vector of system variables and p is a vector of constant system 
parameters, the Jacobian J{x]p) is dehned as the matrix with entries 

T _ dfi{x;p) 

~ dx, • ^ ^ 

Solution to (1) is typically obtained using some iterative method. For 
example, simple Newton’s method approaches solution to (1) by iterating 

Xk+i = Xk - J{xk]p)~^f{xk;p). (3) 

In most of the cases, the model developer needs to provide Jacobian in ad¬ 
dition to model equations. For most engineering problems the governing 
equations are sparsely coupled. That means only a small fraction of Jaco¬ 
bian entries will be nonzero. There is a number of different algorithms for 
solving linear systems that take advantage of the system sparsity to speed 
up computations [7, 8]. To use those algorithms, however, one also needs 
to provide the sparsity pattern. Developing an efficient way for comput¬ 
ing the sparsity pattern of the Jacobian matrix is the key enabling technol¬ 
ogy for solving large-scale nonlinear systems, ordinary differential equations, 
differential-algebraic equations, and optimization problems. Since all of the 
information required to compute Jacobian is contained within model equa¬ 
tions, the computation of the Jacobian and its sparsity pattern can be fully 
automated. 

2.2 Function Derivatives 

There are several ways to compute Jacobian derivatives (2). 

• Compute and implement analytical Jacobian manually. While this will 
lead to fastest numerical computations, it is often not feasible to com¬ 
pute Jacobians manually for large systems. Furthermore, in order to 
compute Jacobian manually, one needs to assume that exact form of 
equations (1) is known a priori. Often, a requirement is to reuse model 
equations in cases where some of the parameters pi are set as variables 
and some of the variables Xi are “hxed” to constant values. Each of 
those cases would require different Jacobians. For large systems, the 
number of all possible combinations would be prohibitive. 

However, if the systems are composed of only a small number of differ¬ 
ent basic components, computing the element stamps for these basic 
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components might be feasible. Circuit simulators, for instance, use this 
approach (cf. [9]). 

• Compute Jacobian numerically. Derivatives (2) can be computed using 
numerical approximation 

^ _ fi{xi,...,Xj + e,...,Xn]p)-fi{x]p) 

Jij ~ ^ 

where e is a small parameter. This approach is general and relatively 
easy to implement, but it is not efficient from computational standpoint 
as it requires additional function evaluations. In addition to that, the 
choice of the approximation parameter £ may affect convergence of the 
solution. If £ is too large, the derivative approximation will be poor. If 
it is too small, too many signihcant digits are lost in the numerator of 
(4) due to the hnite precision of the numerical values. 

• Compute Jacobian symbolically. There is a number of tools and algo¬ 
rithms that can compute derivative expressions given the model equa¬ 
tions (see e.g. [10]). In addition to mathematical algorithms, this 
approach requires equation syntax parsing capability. For large sys¬ 
tems, parsing the equations can be quite time consuming. Parallelizing 
these methods could be quite challenging, as well. 

• Use automatic differentiation to compute Jacobian. In this approach, 
all derivatives are computed automatically at the same time when 
model equations are evaluated. This can be implemented in C-I--I- or 
any other language that supports operator overloading. Jacobian is as¬ 
sembled from available derivatives at the system level. This approach 
does not require any involvement from component model developer. 
Some computational overhead is expected when using automatic dif¬ 
ferentiation. 

In our approach, we selected automatic differentiation as the preferred 
method to compute Jacobian derivatives. Automatic differentiation provides 
exact derivatives without introducing numerical errors and it does not put 
additional burden on component model developers. Since automatic differen¬ 
tiation uses operator overloading, all arithmetic operations and mathematical 
function evaluation will have some overhead. The objective is to have an im¬ 
plementation where such overhead will be small enough so that the overall 
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computational cost is smaller than numerical evaluation of the derivatives. 
Using symbolic manipulations to compute derivatives was ruled out because 
it requires fairly complex equation parsing capability. Furthermore, it is not 
clear how current state-of-the-art symbolic differentiation algorithms scale 
with the size of the problem and if they can be successfully parallelized. 

2.3 Simple Example 

In this subsection, we provide a simple example and illustrate what input 
a sparse nonlinear solver typically requires. Assume we are trying to find a 
steady state solution for a Lorenz system [11]. The residual equations (1) 
can be written as 

a{y-x) = 0 , 

x{p- z)-y = Q, (5) 

xy — /3z = 0. 

Here, a, p, and [3 are constant parameters and x, y, and 2 ; are variables. The 


Jacobian matrix (2) is then given by 


—a a 0 

J = 

p — z —1 —X 


_ y X -I3_ 


Jacobian entry J 13 = 0, while all the other entries have nonzero values, 
generally. The sparsity pattern for Lorenz system is then 

• • 

.... (7) 

• • • 

This tells the solver it does not need to allocate memory for J 13 and perform 
computations with it. Clearly, removing one out of nine Jacobian entries 
does not reduce computational cost signihcantly. Benefits of using sparse 
algorithms will become obvious when we look at larger, real-life problems. 

Analogously, the structure of a system can be represented as a bipartite 
graph, where the bipartite sets of vertices are the equations and variables, 
respectively. Equation fi is then by dehnition connected to variable xj if 
and only if Jij ^ 0. The dependency graph of the Lorenz system is shown 
in Figure 1. These dependency graphs are typically used by Modelica-based 
tools for the causalization of equations and tearing algorithms, see for exam¬ 
ple [ 12 ]. 
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Figure 1: Dependency graph of the Lorenz system. 

3 Sparse Automatic Differentiation 

3.1 Automatic Structure Analysis 

For better clarity, let us hrst discuss sparsity pattern generation alone. Spar¬ 
sity pattern such as the one in (7) is required by the numerical solver at the 
initialization stage to allocate objects required for sparse linear algebra algo¬ 
rithms. During computations, the sparsity pattern (i.e. connectivity struc¬ 
ture) is used by the linear solver to identify structurally nonzero elements 
of Jacobian matrix that enter computation. System connectivity informa¬ 
tion can be used for a number of other analyses such as index reduction 
for differential-algebraic equations (DAEs), partitioning, model causaliza- 
tion, tearing, numerical diagnostics, and many others. These are, however, 
beyond the scope of this paper. 

The approach we propose is to compute residual dependencies on the fly 
along with the residual value computation. To do that, we dehne a math¬ 
ematical object y which is a set containing a real number y and set D, 
which contains integer labels of all dependencies of y. Labels are indepen¬ 
dent variable identihers; typically they are offset values in the solution vector 
as returned by the solver. We denote this object as 

( 8 ) 

For any independent variable x, the corresponding dependency tracking ob¬ 
ject is 

A = {x, (9) 

that is, each independent variable has only trivial self-dependency. Algebraic 
operations on Y, the set of all y, are dehned as follows: 
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• For any C G M and 3^ G Y and algebraic operation * it is 


C*y = {C*y,B>y}. (10) 

• For any two 3^, ^ G Y and mathematical operation * it is 

3 ; * Z = {y * 2 ;, Dj, UDJ. (11) 

• For any function h{y) defined on M, there is a corresponding function 
h{y) defined on Y such that 

h{y) = {h{y), DJ, (12) 

where is the set of dependencies of 3 ^. 

Comparisons between elements of Y are performed with respect to values 
only, disregarding dependencies. For example: 

3^1 > 3^2 yi> y 2 - (13) 

If we define residual equations on Y, rather than M, the residual computation 
will give us both, residual value and the sparsity pattern. Take for example 
Lorenz system (5). The first residual is computed as 

J'l = a ({x, - {y, {rij^}}) 

= - y, {n^} U {ny}} (14) 

= W{x -y),{n^,ny}}, 

and similarly we get 

^2 = {x{p - z) - y,{n^,ny,n^}}, (15) 

^3 = {xy - I3z,{na:,ny,nz}}. (16) 

From dependencies in (14-16), one can obtain sparsity pattern (7) by setting 
Ux = I, Uy = 2, and Uz = 3. Note that this approach for getting sparsity 
pattern is independent of how equations are written. If we, for example, 
write the third residual in (5) as 

/s = M - /3z, (17) 
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where u = xy, then residual evaluation using dependency tracking variables 
gives us 

J 3 = {u, D„} - I3{z, {n^}} 

= {xy,{n^}U{ny}} - {(3z,{n,}} 

= {xy - /3z, {n^, Uy} U {n^}} 

= {xy - /3z, {n^,ny,n^}}, 

which is the same as (16). This property is particularly convenient when 
coding residual equations because it allows reordering computations and us¬ 
ing as many intermediate variables as necessary. Note that all derivatives 
in the equations are uniquely defined in terms of derivatives with respect to 
independent variables. 

3.2 Automatic Differentiation 

To perform sparse automatic differentiation we make a small extension to 
the object we used for sparsity pattern computation. We define Y as a set 
of all 

y = {y, {{n, dny) : n G (19) 

where y and Dy are same as in ( 8 ). Essentially we mapped to each dependency 
a value of partial derivative with respect to that dependency. For independent 
variables 

;^ = {x, {(n,,l)}}. (20) 

Algebraic operations on Y are defined in a similar fashion as in the depen¬ 
dency tracking case: 

• For any C G M, ^ G Y and algebraic operation * defined on Y it is 

C*y = {C*y,{(n, d^(C *y)):ne Dy}}. ( 21 ) 

• For any two y, 2 G Y and algebraic operation * defined on Y it is 

y * £ = {y * z, {(n, dn{y * 2 ;)) : n G Dy U D^}}. (22) 

• For any function h{x) defined on M, there is a corresponding function 
h{y) defined on Y such that 

Ky) = {Ky)^ {{n, h'{y)dny) : n G Dy}}, (23) 

where Dy is the set of dependencies of 3 ^. 
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Comparisons between elements of 3^ are defined in the same way as for de¬ 
pendency tracking data type. 

As an example, let us compute Jacobian derivatives for residual (17). Us¬ 
ing automatic differentiation data type defined in (19), and following algebra 
defined for it, this computation is carried out as 

Js = {m, {(Ua;, dxu), (Uy, dyu)]} - (3 {z, {(n^, 1)}} 

= {xy, y), (uy, a:)}} - {^z, {{n,,(3)]] 

= {xy - /3z, {{n^, y), (uy, x), (n^, -/?)}}. 

Derivatives in make up the third row of the Jacobian matrix (6), when 
Ux = 1, Uy = 2, and = 3. 

4 Prototype Implementation 

Dependency tracking objects like (8) or (19) can be implemented in any pro¬ 
gramming language that supports operator overloading. We created prelim¬ 
inary implementation in C-|--|- mainly for prototyping and testing purposes. 
Here, we outline details of this implementation. 

We create class Variable that stores double precision value and depen¬ 
dency map related to that value. The class overloads all operators dehned 
for the double data type. The prototype implementation is structured like 
this: 

class Variable 
{ 

public : 

Variable (); 

explicit Variable( double value); 

Variable (double value, size_t variablelD); 

Variable( const Variable^ v) ; 

"Variable(); 

// = 

Variable& operator=( const double& rhs); 

Variable& operator=( const Variablefc rhs); 

// += 

Variable& operator+=( const doublefe rhs); 

Variable& operator +=( const Variable& rhs); 

// *= 

Variable& operator *=( const double& rhs); 

Variable& operator *=( const Variablefe rhs); 
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// 


typedef std : : inap<size_t , double> DependencyMap ; 

private : 

double value_; 
size_t variableID_; 
bool isFixed_; 

mutable DependencyMap* dependencies.; 

}; 

This class has several constructors. The constructor that creates Variable 
from double data type is made explicit to prevent possible loss of derivatives 
in accidental implicit data conversions. In addition to variable value and 
variable identifier, the object has boolean flag isFixed_. This flag is used 
when designation of the object needs to change from a variable to a constant 
parameter. The dependency map in this implementation is just a standard 
map between variable identifiers and values of derivatives with respect to 
those variables. Only independent (state) variables have assigned identihers. 
Dependent or temporary variables would obtain their dependency map di¬ 
rectly or indirectly from state variables. For example x, y, and z in equation 
(17) would have identifiers set to 1, 2, and 3, respectively, but u would not 
have an identifier. 

Arithmetic operators are implemented in terms of compound assignment 
operators. For example, here is how *= operator overloading is implemented 
for the Variable class: 

Variablefc Variable:: operator *=( const Variables rhs) 

{ 

// derivation by parts of 0 *this 
scaleDependencies(rhs.value.); 

// compute partial derivatives of rhs and add them to *this 
foreach (DependencyMap::value.type p, *(rhs.dependencies.)) 
(♦dependencies.)[p“>first] += (p->second * value.); 

// compute value of *this 
value. *= rhs.value.; 

return *this ; 

} 

Here, we compute derivative {Ihs ■ rhs)' = {Ihs)' ■ rhs + Ihs ■ {rhs)', where 
Ihs is pointed to by this. Derivatives of Ihs are hrst scaled by the value 
of rhs to obtain the hrst term in the expression for the derivative of the 
product. Then each rhs derivative is multiplied by Ihs value and added to 
corresponding derivatives of Ihs. If the corresponding derivative of Ihs does 
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not exist, a new entry to dependency map is created. The value is computed 
at the end, because it is used in derivative of the product expression. 

The multiplication operator is then implemented as a non-member oper¬ 
ator 

const Variable operator *( const Variablefe Ihs , const Variable& rhs) 

{ 

return Variable(Ihs) *= rhs; 

} 

The copy constructor with Ihs as the argument is used to create the variable 
to be returned and then *= operator is used to multiply that variable by rhs. 

In addition to overloading operators, all mathematical functions from the 
standard C-I--I- library operating on double data type have to be overloaded, 
as well. Take for example sine function whose derivative is (sinr)' = cosx-x': 

namespace std 
{ 

inline Variable sin(const Variablefc x) 

{ 

double val = sin(x.getValue()); 

double der = sin_derivative(x.getValue()); 

Variable res(x); // copy derivatives of x 

res.setValue(val); // set function value f(x) 

res.scaleDependencies(der); // compute derivatives of f(x) 

return res; 

} 

} 

In namespace std we define inline function that takes Variable data type as 
the input. We use copy constructor to retain all derivatives of x, and then 
we multiply them by cosx per chain rule. Value of the sine is computed 
using sine function from the standard library. The function that computes 
the derivative of the sine is defined in the same namespace as Variable type: 

inline double sin_derivative (double x) 

{ 

return std::cos(x); 

} 

The same approach is used for other functions (see Appendix B). 

A simple use of the Variable class is shown in the following code: 

template <typename T> 

void residualFunction(vector<T>& f, 

const vector<T>& x, 
const vector<T>& p) 

const T y = x[0]*x[l]; 

f [0] = p[0]*(x[l] - X [0] ) ; // sigma*(y - x) 

f[l] = x[0]*(p[l] “ x [2] ) - x[l]; // x*(rho - z) - y 
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f [2] = y - p[2]*x[2]; // x*y - beta*z 

} 

int main () 

{ 

const size_t n = 3; 

vector<Variable> x(n), p(n), f(n); 

// decide x, y, and z are variables... 
for (size_t i = 0; i < n; ++i) 

X[i] .setVariableNumber (i) ; 

// ...and sigma, rho, and beta are constant parameters 
for (size_t i = 0; i < n; ++i) 
p[i]. setFixed(true) ; 

// initialize independent variables 
x[0] = 8.0; x[l] = 20.0; x [2] = 2.0/3.0; 

// set constant parameter values 

p[0] = 10.0; p[l] = 8.0/3.0; p[2] = 28.0; 

residualFunction(f, x, p); 
printlncidenceMatrix(f); 
printJacobian(f); 

} 

The function residualFunction computes residual for steady state solution 
of the Lorenz system. Temporary variable y is not really needed in this exam¬ 
ple other than to illustrate how derivative calculation is propagated through 
variables. It is important to note that the Lorenz model in this implementa¬ 
tion does not depend on a specihc data type. Furthermore, the model does 
not assume what are independent variables and what are system parameters. 
This is determined outside the model. In this example, elements of vector 
X are set to be independent variables and elements of vector p constant pa¬ 
rameters in the main function. This is consistent with the problem dehned 
in (5) where we look for a steady state solution given parameters a, p and (3. 

One residual evaluation with Variable data type also computes sparsity 
pattern and Jacobian derivatives. Function printIncidenceMatrix outputs 
sparsity pattern: 

A = [1 1 0; 

1 1 1; 

111]; 

Function print Jacobian outputs: 


J = [-10 10 0; 

2 -1 - 8 ; 


14 


20 8 -28]; 


Outputting the matrices in this format enables postprocessing in Matlab to 
detect, for instance, singular or ill-conditioned Jacobians and to easily plot 
results. For simplicity, we omit implementation of the two output functions. 

It is also possible to consider elements of the vector p as system variables 
and elements of vector x as constant system parameters. All one has to do 
is to set: 

// decide x, y, and z are constant parameters... 
for (size_t i = 0; i < n; ++i) 
x[i] .setFixed(true); 

// ...and sigma, rho, and beta are variables 
for (size_t i = 0; i < n; ++i) 
p[i] .setVariableNumber (i) ; 

In this problem we look for parameters a, p, and /3 such that a hxed point 
solution of (5) is at x = 8, p = 20, and z = 2/3. System sparsity pattern 
and Jacobian in that case are obtained as 

A = [100; 

0 10 ; 

0 0 1 ]; 


and 


J = [12 0 0; 

0 8 0 ; 

0 0 -0.667]; 

respectively. This change can be made at runtime without need to recompile 
the system model. Any other selection of system parameters and variables 
can be made in the same way. Vectors x and p are used merely to denote 
nominal system variables and parameters. 

Typically, one would run residual evaluation during solver initialization to 
get sparsity pattern, and then run residual evaluation every time Jacobian 
is required during solver iterations. Note that residual vector in this case 
is in fact an implementation of a compressed row sparse matrix. The only 
difference is that each row, in addition to matrix elements, also holds a 
corresponding residual vector element. 
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5 Preliminary Benchmarking Results 

Dependency tracking and automatic differentiation algorithms as implemented 
here are exact. The main challenge is to ensure that computational overhead 
introduced by automatic differentiation scales well with the size of the sys¬ 
tem. To provide a preliminary assessment of how computational cost for 
automatic differentiation scales with the size of the model (i.e. number of 
model equations), we perform several numerical experiments. 



Figure 2: Benchmarking test case: A simple electrical grid model. 

As a benchmark problem, we select a simple microgrid model as shown 
in Figure 2. A single alternating current (AC) generator is connected to a 
rectiher that converts power to direct current (DC) and supplies it to a DC 
bus. Several passive AC loads are connected to the DC bus, each through 
a separate inverter that convert DC power from the bus to 60 Hz AC power 
required by the load. The size of this system can be easily scaled up by 
simply adding more loads to the bus. System parameters are set so that 
simulation results are “self-validating” - voltage at each load is the same 
as the voltage produced by the generator: 100 V, 60 Hz sinusoidal. A more 
detailed description of the test case is given in Appendix A. 

The electrical grid model is cast in form of differential-algebraic equations. 
We simulate hrst 0.1 s of the grid operation using the Rythmos package from 
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the Trilinos library. We use an implicit variable-order variable-stepsize back¬ 
ward differentiation method to solve the equations [13]. Nonlinear solution 
at each time step is obtained using a sparse direct method. The method 
requires residual equations and system Jacobian to be provided. We mea¬ 
sure overall computation time and model evaluation time. In our case, model 
evaluation is a residual evaluation using Variable data type. Jacobian is eval¬ 
uated automatically together with residual at each model evaluation call per 
design of Variable class. 

For differential-algebraic equations, residual equations are of the form 

F{v,v,t) = 0. (24) 


Unknown variables v in our case are typically node voltages. Jacobian for 
differential-algebraic equations is somewhat different than the one defined in 
(2). It is typically given as 


^ dF dF 

J — , 

ov ov 


(25) 


where parameter a is provided by the solver and is related to numerical 
integration scheme used in simulation. Example of the Jacobian sparsity 
pattern for the benchmarking test case is shown in Figure 3. 



Figure 3: Sparsity pattern for the grid model with 30 loads. Only 1.8% of 
Jacobian elements are structurally nonzero. 

In the serial case, we ran simulations for grids with 100-600 loads. For 
these simulations it takes roughly 8,000 integrator steps and 50,000-60,000 
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function evaluations to complete regardless of the size of the system. We 
hnd that average computational time per call grows linearly with the size 
of the system as shown in Figure 4. Data points for the computational cost 
of function evaluation £t particularly well to the linear £t. Furthermore, 
function evaluation becomes smaller fraction of the overall cost as the size of 
the system increases. 



Figure 4: Computational time per call for integration step and function 
evaluation. 

We compare average cost of function evaluation for our sparse automatic 
differentiation prototype with dense automatic differentiation using Sacado 
package from the Trilinos library. The prototype uses map from standard 
C++ library (O(logn) cost), whereas Sacado uses dense vector (0(1) cost) to 
store and access derivatives. Sparsity pattern information is provided to ex¬ 
ternally in Sacado case, so that only structurally nonzero derivatives are com¬ 
puted. For systems of this size, it was expected that dense algorithm would 
outperform the sparse automatic differentiation. Both algorithms evaluate 
the same derivatives and the dense approach has faster access to derivatives. 
The only downside of the dense approach is that it has to allocate larger 
chunks of memory to store derivatives. Yet, our results suggest that memory 
management alone may cause computational cost to grow quadratically with 
the size of the system when using dense automatic differentiation (Figure 5). 

Parallelizing simulations that use sparse automatic differentiation to com¬ 
pute Jacobian is fairly straightforward in MPI framework when using our 
approach. For testing purposes, we implemented a simple parallelization 
scheme where generator and rectifier are simulated on one node and simu¬ 
lations of inverters and loads are evenly distributed over remaining nodes. 
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Figure 5: Function evaluation cost when using sparse and dense automatic 
differentiation. 

We show here results of an MPI simulation on 16 CPU cores (4 nodes). The 
number of loads in the system was varied from 1,000-4,000 (roughly 10,000- 
30,000 equations). The results show again linear scaling as the size of the 
system increases as shown in Figure 6. Same as in the serial case, the cost of 
the function evaluation grows slower than the cost of the solver as the size 
of the system increases. 



Figure 6: Computational time per call for integration step and function 
evaluation. 
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6 Conclusion 


Our analysis and benchmarking results suggest that using abstract elemen¬ 
tary algebra approach for sparse automatic differentiation is a promising 
direction. Linear scaling of the computational cost and the ease of paral¬ 
lelization indicates that this approach is particularly suitable for massively 
parallel computations. Overhead of using abstract elementary algebra be¬ 
comes a small fraction of the overall computational time for large systems. 

The prototype implementation leaves room for code optimization. The 
residual vector implemented in terms of Variable data type is de facto a 
compressed row sparse Jacobian matrix with the extension that each row is 
associated with corresponding residual value. Intermediate variables, such as 
u in (17), can be understood as sparse matrix rows, as well. However, they are 
not part of the Jacobian, they are used just to complete the chain rule. Typ¬ 
ically, intermediate variables are used by system modelers for convenience 
to write model equations in a more compact form. In the current imple¬ 
mentation, derivatives of residual functions are reallocated at every function 
evaluation based on dependency tracking mechanism. Since this information 
does not change during solver iterations, dependency structure of the system 
could be precomputed once and then reused at subsequent solver iterations. 
This could be done easily for residual vectors, which are typically passed by 
reference to models. Intermediate variables are typically local variables used 
by the modeler to simplify the equations, so they may have to be reallocated 
at every iteration anyways. 

Temporary scope of residual vector and Jacobian elements might be suit¬ 
able for matrix-free methods [14], where one wants to avoid storing the entire 
matrix. 

When simulating our test cases, we had to copy residual from the vec¬ 
tor of variables to the Epetra vector and Jacobian derivatives to Epetra 
compressed-row sparse matrix, per requirements of the Rythmos solver. This 
added small additional overhead to the computation. The full power of the 
proposed approach could be demonstrated with numerical solvers that do 
not require specihc data formats, but instead provide abstract interfaces to 
all linear and elementary algebra operations. While such solvers are still not 
part of the mainstream, a lot of activities have been done in that direction 
as for example in Tpetra project [6]. 

Using abstract elementary algebra has potential applications way beyond 
automatic differentiation. As we have shown in this paper, it could be used 
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for reconfiguring models at runtime. Constant parameters could be changed 
into variables and vice versa. Abstract data types also could be used for 
diagnostics, for example to identify structurally singular Jacobians. Fur¬ 
thermore, this approach could be used for preprocessing model equations for 
index reduction of differential-algebraic equations or tearing algorithms for 
system decomposition. 

Abstract elementary algebra can also help code reuse. Same model code 
can be reused for local sensitivity analysis or embedded uncertainty quantih- 
cation, simply by using different template parameter. The same holds true for 
solvers that provide abstract interfaces for elementary algebra. More reuse 
streamlines code verihcation and improves development efficiency, which are 
critical for any large scale computation. All of these will be pursued in sub¬ 
sequent work. 
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A Test Case Description 

Electrical grid in Figure 2 is a common motif in power systems. Since math¬ 
ematical modeling methods for power systems are not commonly known out¬ 
side electrical engineering community, we provide here brief overview of the 
governing equations used in our test case. For a detailed description of the 
component models, the reader should consult for example [15] and refer¬ 
ences therein. The equations are derived and the system is composed using 
modihed nodal analysis approach [16]. 

We assume the generator shown on the left side in Figure 2 is an ideal 3- 
phase generator and set residual equations for voltages on generator terminals 
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a, b, and c as: 

(26) 

(27) 

(28) 


0 = Vga- V()Sin(i:nt), 

0 = Vgb — Vosin(a;t + 27r/3), 
^ = Vgc — Vq sin(a;t + dvr/S). 



Figure 7: Schematic of the rectiher model. 

The generator is connected to a rectiher, which converts 3-phase AC 
power to DC power. Rectiher and hlter schematics are shown in the Figure 
7. Kirchhoh’s current law for the rectiher can be cast in terms of residual 
equations as 


0 ^ga ^ni^ga ^/) T ^ga)y (^^) 

0 igb loi^gb ^/) T ^oiVn '^gb)i (^0) 

0 igc loiVgc ^f') T ^niVn Vgc)- (31) 


Here, iga, igb, and igc are generator phase currents entering rectiher nodes a, 
b, and c; u/ is the voltage at node / (rectiher’s positive terminal, connected 
to the hlter) and Vn is the node voltage at rectiher’s negative terminal. 
Current through diodes Id{v) in the rectiher is modeled as 


Id{v) 




In our simulations we chose temperature T = 300 K, ideality factor n = 2, 
and saturation current to be A = 18.8 nA. Boltzmann constant k ^ 1.38 x 
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10 J/K. Voltage v in the diode current function is the difference between 


node voltages at diode’s anode and cathode terminals. 

Residual equations for the hlter can be written as 

0 = loivga - Vf) + loivgb - Vf) + Inivgc “ Vf) - ii, (32) 

0 = (j)L- Lii, (33) 

0 = 0L - (u/-Up), (34) 

Q = qc - C{vp-Vn), (35) 


where il and (j)L are inductor’s current and flux, respectively, L is the in¬ 
ductance, qc is charge on the capacitor, C is the capacitance, and Vp is the 
voltage at positive terminal of the DC bus. 

Kirchhoff’s current law gives following equations for the DC bus (Fig¬ 
ure 2): 

0 = - qc + ib 

k 

0 = loiVn - Vga) + IniVn “ Vgb) + loiVn “ Upc) 
k 

where and are a-phase of the voltage across and current through the 
load k {a = a,b,c). Functions d^\t) are inverter modulation signals that 
describe 3-phase AC waveform at inverter outlet. We choose modulation 
signals to be the same for all inverters and produce sinusoidal output: 

— ?7isin(a;t), 

= ''^sin(a;f -|- 27r/3), 

= msin(a;f -|- 47r/3). 

Here we set u to be the same as the frequency of the generator. Furthermore, 
we set 

27r 


so that voltage amplitude at each load is the same as the generator voltage 
amplitude. This choice was made merely for model verihcation convenience 
(again, interested readers are referred to reference [15] for more details). 


(36) 

(37) 
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Ua 

hh 

Uc 


Figure 8 : Schematic of the inverter model with averaged pulse width modu¬ 
lation. 


The inverter model we used is an averaged model without details of pulse 
width modulation. The equivalent circuit model of such inverter model is 
shown in Figure 8 . The model consists of three ideal current-controlled cur¬ 
rent sources and three ideal voltage-controlled voltage sources (one for each 
phase). The current sources are controlled by the load currents as 

Zq, ® 0 , 6 , c, 

where Pq, is load current (Figure 8 ). Voltage sources are controlled by the 
DC bus voltage as 


{t) {Vp-Vn), a = tt, 6, C. 

Here, Vp and are node voltages on positive and negative terminals of the 
DC bus, respectively. By using Kirchhoff’s current law, the equations for the 
inverter are obtained as 


0 = 4? - ' 

0 = - G^{vP 

0 = 4 ? - 

0 = - t.f 


- -.i') + 

- fo*') + dP{t)Gt(vlf> - tij,*’), 


) + Gkiv^ * — + Gk{v 


(ty. 


c 


.m, 


(38) 

(39) 

(40) 

(41) 


Using Kirchhoff’s voltage law, we obtain equations for voltages at load ter- 
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minals Via- 


(k) 1 + d^\t) . . . . 

0 = vi ^ -Vn - {Vp -Vn), (42) 

0 = 4^^ -Vn - {Vp - Vn), (43) 

(k) 1 + dlcHt) , s , , ,s 

0 = v\l^ -Vn - 

The entire system is described by residual equations (26-44). There are 
12 + 7A^ system variables, where N is the number of AC loads connected 
to the bus. System variables are generator node voltages Vga, Vgb, and Vgc, 
generator currents iga, igb, and igc] rectiher node voltage u/; DC bus positive 
and negative voltages Vp and hlter internal variables - inductor current 
ijb, inductor flux (f)^, and capacitor charge qc] currents through ideal voltage 
sources in inverter model , and and load node voltages Vi^\ v\l\ 

vf^, and Index k = 1,...,A^ denotes a load. In our tests, we set 

parameters to the following values: Generator frequency u = 27r60 rad/s, 
load conductances are all equal and set to G = 0.01 S, capacitance in the 
hlter is G = 0.1 mF, and inductance in the hlter is L = 20 mH. 


B Variable Class Code 

B.l Class Declarations and Inline Implementation 


/* ! 

\file Variable.hpp 

*/ 

#ifndef SAD_VARIABLE_HPP 
#define SAD_VARIABLE_HPP 

#include <map> 

#include <vector> 

#include <cstdlib> 

#include <cmath> 

#include <fstreain> 

#include <cassert> 

#include <liinits> 

#include <iostreain> 

#include <iOmanip > 

#include <string> 

#include <boost/foreach.hpp> 
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#include <boost/math/spe cial_functions/sign.hpp > 


namespace SAD 

{ 


\brief The Variable class is used to store the unknowns of the 
system and define abstract elementary algebra. 

\author Stefan Klus 
\author Slaven Peles 

*/ 

class Variable 

{ 

public : 

/ * ! 

\brief Default constructor. 

*/ 

Variable () 

: value_ (0.0) , 

variableNumber_(INVALID_VAR_NUMBER) , 
isFixed_ (false) , 

dependencies_ (new DependencyMap) 

{ 

} 

/ * ! 

\brief Constructor which initializes the value. 

*/ 

explicit Variable( double value) 

: value_(value), 

variableNumber_(INVALID_VAR_NUMBER) , 
isFixed_ (false) , 

dependencies_ (new DependencyMap) 

> 

/ * ! 

\brief Constructor which initializes the value and variable 
number. 

*/ 

Variable( double value, size_t variableNumber) 

: value_(value), 

variableNumber_(variableNumber), 
isFixed_ (false) , 

dependencies_ (new DependencyMap) 

{ 

(*dependencies_)[variableNumber_] = 1.0; 

} 

/ * ! 

\brief Copy constructor. 

*/ 

Variable( const Variable& v) 

: value_(v.value_), 


26 


variableNumber.(INVALID_VAR_NUMBER) , 
isFixed_ (false) , 

depeiidencies_(new DependeiicyMap(*v.dependencies_)) 

{ 

} 

/ * ! 

\brief Destructor deletes the dependency map. 

*/ 

"Variable() 

delete dependencies_; 

} 

/ * ! 

\brief Assignment operator. Assigning double value to 
Variable removes its dependencies. Use only if you know 
what you are doing. 

*/ 

Variable^ operator=( const doubled rhs) 

{ 

value_ = rhs; 

dependencies_->clear(); 
return * this ; 

} 

/ * ! 

\brief Assignment operator. 

This operator: 

- assigns value from the right hand side 

- leaves variable ID unchanged 

- clears any existing and adds new dependencies from rhs 

*/ 

Variable& operator=( const Variable& rhs) 

{ 

if (this == &rhs) // self-assignment 
return *this ; 

// set value from rhs 
value_ = rhs.value_; 

// if rhs is a constant so is *this 
setFixed(rhs . isFixedO) ; 

// set dependencies from rhs 

dependencies_->clear (); // clear map just in case 

addDependencies(rhs); // use only dependencies from the rhs 

return * this ; 

} 


/ * ! 

\brief Operator () returns the value of a variable. 
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This is just short notation to avoid using 
getValue and setValue. 

*/ 

doublefe operatorOO 
return value_; 

} 

/ * ! 

\brief OperatorO returns the value of a variable 
(const version). 


This is just short notation to avoid using getValue. 

*/ 

const doublefc operatorOO const 


return value_; 

} 


/ * ! 

\brief Return the current value of the variable. 

*/ 

double getValueO const 
return value_; 

} 


/ * ! 

\brief Overwrite the current value of the variable. 

*/ 

void setValue( double value) 

{ 

value_ = value; 

} 


/ * ! 

\brief Return derivative of *this with respect to 
dependency i. 

*/ 

const double der(size_t i) const 
return (*dependencies.)[i]; 

} 


/ * ! 

\brief Returns the variable number. 

This number is assigned to state variables (variables 
updated directly by the solver) only. 

*/ 

size.t getVariableNumber () const 
return variableNumber.; 


\brief Sets the variable number. 

*/ 

void setVariableNumber(size_t variableNumber) 

dependencies_->clear(); 

variableNumber_ = variableNumber; 

(*dependencies_)[variableNumber_] = 1.0; 

} 

/ * ! 

\brief Checks whether the variable was registered as 
an unknown of the system. 

INVALID_VAR_NUMBER is used to mark parameters and 
temporary variables 

*/ 

bool isRegistered () const 

return variableNumber. != INVALID_VAR_NUMBER; 

} 

/ * ! 

\brief Checks whether the variable is fixed or not. 

*/ 

bool isFixedO const 

{ 

return isFixed_; 

} 


/* 


*/ 


\brief Turns variable into parameter, or vice versa. 


void setFixed (bool b 

{ 

isFixed. = b; 

} 


false ) 


// get the ’ input set ’ of a variable 

typedef std : :map<size_t , double> DependencyMap; 

inline const DependencyMap& getDependencies() const; 

// set as the independent state variable and assign ID to it 
inline void registerVariable(std::vector<Variable*>& x, 

const size_t& offset); 

// adds all dependencies of v to *this 

inline void addDependencies( const Variablefe v); 

// scale dependencies (derivatives) by scalar Oa c. 
inline void scaleDependencies (double c); 

// print to output stream 

inline void print(std::ostreamfe os) const; 

// += 

inline Variable& operator+=( const double& rhs); 
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inline Variable& operator + =( const Variablefe rhs) : 

// -= 

inline Variablefe operator -=( const double& rhs); 
inline Variable& operator -=( const Variablefe rhs): 

II * = 

inline Variable& operator *=( const double& rhs); 
inline Variable& operator *=( const Variablefe rhs): 


// / = 

inline Variable& operator /=( const double& rhs); 
inline Variable& operator /=( const Variablefe rhs): 


private : 

double value.; ///< Value of the variable, 

size.t variableNumber.; ///< Independent variable ID 

bool isFixed.; ///< Constant parameter flag. 

mutable DependencyMap* dependencies.; 
static const size.t INVALID.VAR.NUMBER = -1; 

}; 


II - 

// non-member operators and functions 
//- 


// unary - 
inline const 

Variable 

operator “(const 

// + 

inline 

const 

Variable 

operator+( const 

inline 

const 

Variable 

operator+( const 

inline 

const 

Variable 

operator+( const 

// - 

inline 

const 

Variable 

operator“(const 

inline 

const 

Variable 

operator “(const 

inline 

const 

Variable 

operator “(const 

// * 

inline 

const 

Variable 

operator*(const 

inline 

const 

Variable 

operator*(const 

inline 

const 

Variable 

operator*(const 

// / 

inline 

const 

Variable 

operator/(const 

inline 

const 

Variable 

operator/(const 

inline 

const 

Variable 

operator/(const 


Variablefc v); 

Variablefe Ihs, const Variablefe rhs); 
Variablefe Ihs, const double& rhs); 
double& Ihs, const Variable& rhs); 

Variablefc Ihs, const Variable& rhs); 
Variablefe Ihs, const doubled rhs); 
double& Ihs, const Variable^ rhs); 

Variablefe Ihs, const Variablefe rhs); 
Variablefe Ihs, const doubled rhs); 
double& Ihs, const Variable^: rhs); 


Variablefe Ihs, const Variable& rhs); 
Variablefe Ihs, const double& rhs); 
double& Ihs, const Variable^ rhs); 


// == 
inline 
inline 
inline 


const bool 
const bool 
const bool 


operator ==( const Variablefe Ihs, const Variable^ rhs); 
operator ==( const Variable& Ihs, const doublefe rhs); 
operator ==( const doubled Ihs, const Variable& rhs); 
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// ! = 


inline 

const 

bool 

operator !=(const 

Variable& 

Ihs 

, const Variablefe rhs) 

inline 

const 

bool 

operator !=(const 

Variablefe 

Ihs 

, const double& rhs); 

inline 

const 

bool 

operator !=(const 

doubled Ihs, 

const Variable& rhs); 

// < 
inline 

const 

bool 

operator <( const 

Variable& 

Ihs , 

const Variablefe rhs); 

inline 

const 

bool 

operator <( const 

Variable& 

Ihs , 

const double& rhs); 

inline 

const 

bool 

operator <( const 

doublefe Ihs, const Variable& rhs); 

// > 
inline 

const 

bool 

operator >( const 

Variable& 

Ihs , 

const Variablefe rhs); 

inline 

const 

bool 

operator >( const 

Variable& 

Ihs , 

const double& rhs); 

inline 

const 

bool 

operator >( const 

doublefc Ihs, const Variable& rhs); 

// <= 
inline 

const 

bool 

operator <=( const 

Variablefe 

Ihs 

, const Variable^ rhs) 

inline 

const 

bool 

operator <=( const 

Variablefe 

Ihs 

, const double& rhs); 

inline 

const 

bool 

operator <=( const 

doubled Ihs, 

const Variablefc rhs); 

// >= 
inline 

const 

bool 

operator >=( const 

Variable& 

Ihs 

, const Variable& rhs) 

inline 

const 

bool 

operator >=( const 

Variable& 

Ihs 

, const doublefe rhs); 

inline 

const 

bool 

operator >=( const 

doubled Ihs, 

const Variable& rhs); 


inline std::ostreamfe operator <<(std::ostream& os, const Variable^ v); 
inline Variable& operator <<(Variable& u, const Variable& v); 

inline std::istreamfe operator >>(std::istream& is, Variable& v); 


y II namespace SAD 

#include "VariableImplement at ion.hpp" 

#include "VariableOperators.hpp" 

#endif // SAD_VARIABLE_HPP 

B.2 Implementation of Member Operators 


/* ! 

\file Variable Implement at ion.hpp 

*/ 

#ifndef SAD_VARIABLE.IMPLEMENTATI0N_HPP 
#define SAD.VARIABLE.IMPLEMENTATION.HPP 

#define foreach BOOST.FOREACH 


namespace SAD 

{ 

/* ! 

\brief Returns the list of derivatives. 

*/ 

const Variable::DependencyMapi: Variable::getDependencies () const 
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{ 

assert(dependencies. != 0); 
return ^dependencies.; 

} 


/ * ! 

\brief Registers a variable as an unknown of the system and 
adds a pointer to the global @a x vector. 

*/ 

void Variable :: registerVar iable ( std :: vector < Variable *>& x, 

const size.t& offset) 

{ 

setVariableNumber(offset); // define global variable number 

setFixed( false ); // not a constant 

X [offset] = this ; 

} 


/ * ! 

\brief Adds all dependencies of v to *this. 

*/ 

void Variable::addDependencies( const Variablefc v) 

{ 

foreach (DependencyMap::value.type p, *(v.dependencies.)) 
(*dependencies.)[p.first] = p.second; 

} 


/* ! 

\brief Multiplies each partial derivative of @a this by @a c. 

*/ 

void Variable::scaleDependencies (double c) 

{ 

foreach (DependencyMap::value.type p, ^dependencies.) 

(*dependencies.)[p.first] *= c; 

} 


/* ! 

\brief Prints the value and input set of the variable. 

*/ 

void Variable::print(std::ostreamfe os) const 

{ 

os << value.; 

if (isFixed.) 

{ 

os << "u(fixed)"; 
return ; 

} 

if (variableNumber. != INVALID.VAR.NUMBER) 

{ 

os << "u(variableu" << variableNumber. << ")"; 
return ; 


32 


> 


if (dependeiicies_ ! = NULL && ! depeiidencies_->empty ()) 

os << "udependencies: u Cu" ; 

foreach (DependencyMap::value_type p, ^dependencies.) 

os << "(" << p.first << ",u" << p.second << ")u": 

os << ; 

} 

} 


H - 

// compound assignment operators 
//- 

/ * ! 

\brief Compound addition-assignment operator. Right hand 
side is a built-in double type. 

*/ 

Variable& Variable:: operator+=( const double& rhs) 

{ 

value. += rhs; 
return * this ; 

} 

/* ! 

\brief Compound addition-assignment operator. Right hand side 
is Variable type. 

*/ 

Variable& Variable:: operator+=( const Variablefe rhs) 

{ 

// compute partial derivatives of *this 

foreach (DependencyMap::value.type p, *(rhs.dependenties.)) 

(*dependenties.)[p.first] += (p.second); 

// compute value of *this 
value. += rhs.value.; 

return *this ; 

} 

// -= 

/* ! 

\brief Compound substruction-assignment operator. Right hand 
side is a built-in double type. 

*/ 

Variable& Variable:: operator -=( const double& rhs) 

{ 

value. -= rhs; 
return *this ; 

} 

/ * ! 

\brief Compound substruction-assignment operator. Right hand 
side is a Variable type. 

*/ 
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Variable& Variable:: operator -=( const Variablefe rhs) 

{ 

// compute partial derivatives of *this 

foreach (DependencyMap::value_type p, *(rhs.dependencies_)) 

(*dependencies_)[p.first] -= (p.second); 

// compute value of *this 
value_ -= rhs.value_; 

return * this ; 

} 

// *= 

/ * ! 

\brief Compound multiplication - assignment operator. Right hand 
side is a built-in double type. 

*/ 

Variable& Variable:: operator *=( const doublefe rhs) 

{ 

// Compute derivatives of this 
scaleDependencies(rhs); 

// compute value 
value_ *= rhs; 

return *this ; 

} 

/ * ! 

\brief Compound multiplication - assignment operator. Right 
hand side is a Variable type. 

*/ 

Variable& Variable:: operator *=( const Variablefe rhs) 

{ 

// derivation by parts of 0 *this 
scaleDependencies(rhs.value_); 

// compute partial derivatives of rhs and add them to *this 
foreach (DependencyMap::value_type p, *(rhs.dependencies_)) 

(*dependencies_)[p.first] += (p.second * value_); 

// compute value of this 
value_ *= rhs.value_; 

return *this ; 

} 

// /= 

/ * ! 

\brief Compound division-assignment operator. Right hand side 
is a built-in double type. 

*/ 

Variable& Variable:: operator /=( const doublefe rhs) 

{ 

double inverseRhs = 1.0/rhs; 

// compute derivatives of Oa *this 


34 



scaleD 

ependencies(inve 

rseRhs 

) ; 




value. 

*= inverseRhs; 






return 

* this ; 





} 







/ * ! 








\brief 

Compound divisi' 

on - as s 

ignment 

operate 

r. Right hand 


side i 

s a Variable typ 

e . 




*/ 







Var 

iable& 

Variable:: operat 

or /=( c 

onst Va 

riablefc 

rhs) 

{ 








double 

inverseRhs 

= 1 

.0/rhs. 

value.; 



double 

inverseRhsSq = 

invers 

eRhs * 

inverseRhs; 


// der 

ivation by parts 

of @a 

♦ this 




scaleD 

ependencies(inve 

rseRhs 

) ; 




foreach (DependencyMap 

::valu 

e type 

p, ♦(rhs 

.dependencies.)) 


(* 

dependencies.)[p 

.first 

] -= (p 

.second 

* value. ♦ inverseRhsSq) 


// com 

pute value of th 

is 





value. 

*= inverseRhs; 






return 

* this ; 






} 

} // namespace SAD 

#endif // SAD_VARIABLE_IMPLEMENTATION_HPP 


B.3 


Implementation of Non-member Operators 


/ * ! 

\file VariableOperators.hpp 

*/ 

#ifndef SAD.VARIABLE.OPERATORS_HPP 
#define SAD.VARIABLE.OPERATORS.HPP 


// Define derivatives of standard library mathematical 
// These definitions could be used to add smoothing in 
// derivatives are discontinuous, 
namespace SAD 
{ 

// - 

// non-member operators and functions 
11 - 


functions here, 
cases where 


// unary - 

const Variable operator -( const Variablefc v) 

{ 

return -1.0* v; 

} 


// + 


35 




const Variable operat or + ( const Variablefe Ihs, const Variablefe rhs) 

{ 

return Variable(Ihs) += rhs; 

} 

const Variable operat or + ( const Variablefe Ihs, const doublei: rhs) 

{ 

return Variable(Ihs) += rhs; 

} 

const Variable operat or + ( const doublefe Ihs, const Variableft rhs) 

{ 

return Variable(rhs) += Ihs; 

} 

// - 

const Variable operator -( const Variablefe Ihs, const Variablefe rhs) 

{ 

return Variable(Ihs) -= rhs; 

} 

const Variable operator -( const VariableSt Ihs, const doubleft rhs) 

{ 

return Variable(Ihs) -= rhs; 

} 

const Variable operator -( const doublefe Ihs, const Variablei: rhs) 

{ 

return Variable(Ihs) -= rhs; 

} 

// * 

const Variable operat or *( const Variablefe Ihs, const Variablefe rhs) 

{ 

return Variable(Ihs) *= rhs; 

} 

const Variable operat or *( const Variablefe Ihs, const doubleS: rhs) 

{ 

return Variable(Ihs) *= rhs; 

} 

const Variable operat or *( const doublefe Ihs, const Variableft rhs) 

{ 

return Variable(Ihs) *= rhs; 

} 

// / 

const Variable oper at or /( const VariableSt Ihs, const VariableSt rhs) 

{ 

return Variable(Ihs) /= rhs; 

} 

const Variable oper at or /( const VariableSt Ihs, const doubleSt rhs) 

{ 

return Variable(Ihs) /= rhs; 

} 
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const Variable operator /( const double& Ihs, const Variable^ rhs) 


return 

} 

Variable(Ihs) /= rhs; 

// == 

const bool 

{ 

return 

} 

operator ==( const Variable& Ihs, const Variablefc rhs) 

Ihs.getValue() == rhs.getValue(); 

const bool 

{ 

return 

} 

operator ==( const Variable^ Ihs, const double& rhs) 

Ihs.getValue () == rhs; 

const bool 

{ 

return 

} 

operator ==( const doublefc Ihs, const Variable& rhs) 

Ihs == rhs.getValue(); 

// ! = 

const bool 

{ 

return 

} 

operator !=( const Variables Ihs, const Variablefc rhs) 

!(Ihs.getValue() == rhs.getValue()); 

const bool 

{ 

return 

} 

operator !=( const Variable& Ihs, const double& rhs) 

! (Ihs.getValue ( ) == rhs); 

const bool 

{ 

return 

} 

operator !=( const doublefc Ihs, const Variable& rhs) 

!(lhs == rhs.getValue()); 

// < 

const bool 

{ 

return 

} 

operator <( const Variable& Ihs, const Variable& rhs) 

Ihs.getValue() < rhs.getValue(); 

const bool 

{ 

return 

} 

operator <( const Variable^ Ihs, const double& rhs) 

Ihs.getValue() < rhs; 

const bool 

{ 

return 

} 

operator <( const doublefc Ihs, const Variable& rhs) 

Ihs < rhs.getValue(); 

// > 

const bool 

operator >( const Variable^ Ihs, const Variable& rhs) 


{ 
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return 

} 

Ihs.getValue() > rhs.getValue(); 

const bool 

{ 

operator >( const Variable^ Ihs, const double& rhs) 

return 

} 

Ihs.getValue() > rhs; 

const bool 

{ 

operator >( const doublefe Ihs, const Variable& rhs) 

return 

} 

Ihs > rhs.getValue(); 

// < = 

const bool 

{ 

return 

} 

operator <=( const Variable& Ihs, const Variablefc rhs) 

Ihs < rhs 1 1 Ihs == rhs ; 

const bool 

{ 

operator <=( const Variable^ Ihs, const double& rhs) 

return 

} 

Ihs < rhs 1 1 Ihs == rhs ; 

const bool 

{ 

operator <=( const doublefe Ihs, const Variablefe rhs) 

return 

} 

Ihs < rhs 1 1 Ihs == rhs ; 

// > = 

const bool 

{ 

return 

} 

operator >=( const Variable^ Ihs, const Variablefe rhs) 

Ihs > rhs 1 1 Ihs == rhs ; 

const bool 

{ 

operator >=( const Variable& Ihs, const double& rhs) 

return 

} 

Ihs > rhs 1 1 Ihs == rhs ; 

const bool 

{ 

operator >=( const doublefc Ihs, const Variable& rhs) 

return 

} 

Ihs > rhs 1 1 Ihs == rhs ; 

//- 



// non-member operators 


//- 


/* 

\brief 

*/ 

Stream insertion operator for variables. 


std::ostreamA: operator <<(std::ostream& os, const Variablefc v) 

{ 

V.print(os); 
return os; 

} 
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/* 

\brief Adds all dependencies of a variable, i.e. v 
that V depends on u. 

\attention Use only if you know what you are doing! 

*/ 

Variable& operator <<(Variable& u, const Variable& v) 

{ 

u.addDependencies(v) ; 
return u; 

} 


/* 


\brief Stream extraction operator for variables. 

*/ 

std::istreamA: operator >>(std::istream& is, Variablefc 

{ 

return is >> v(); 


} 


// 

// 

// 


definitions 


of 


cmath functions derivatives 


V 


III Derivative of sine 

inline double sin_derivative (double x) 

{ 

return std::cos(x); 

} 

III Derivative of cosine 

inline double cos_derivative (double x) 

return -std::sin(x); 

} 

III Derivative of tangent 

inline double tan_derivative (double x) 

{ 

return 1.0/(std::cos(x)*std::cos(x)); 

} 

III Derivative of arc sine 

inline double asin_derivative (double x) 

{ 

return std::sqrt(1.0 - x*x); 

} 

III Derivative of arc cosine 

inline double acos_derivative (double x) 

{ 

return -std :: sqrt(1.0 - x*x); 

} 


<< u means 
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Ill Derivative of arc tangent 

inline double atan_derivative (double x) 

{ 

return 1.0/(1.0 + x*x); 

} 

III Derivative of hyperbolic sine 
inline double sinh_derivative (double x) 

{ 

return std::cosh(x); 

} 

III Derivative of hyperbolic cosine 
inline double cosh_derivative (double x) 

{ 

return std::sinh(x); 

} 

III Derivative of hyperbolic tangent 
inline double tanh_derivative (double x) 

{ 

return 1.0/(std::cosh(x)*std::cosh(x)); 

} 

III Derivative of hyperbolic arc sine 
inline double asinh_derivative (double x) 

{ 

return std::sqrt(1.0 + x*x); 

} 

III Derivative of hyperbolic arc cosine 
inline double acosh_derivative (double x) 

{ 

return std::sqrt(x*x - 1.0); 

} 

III Derivative of hyperbolic arc tangent 
inline double atanh_derivative (double x) 

{ 

return 1.0/(1.0 - x*x); 

} 

III Derivative of exponential function, 
inline double exp_derivative (double x) 

{ 

return std::exp(x); 

} 

III Derivative of natural logarithm function, 
inline double log_derivative (double x) 

{ 

return 1.0/x; 

} 

III Derivative of logarithm to base 10 function: 1/(x*log (10)). 
inline double loglO_derivative (double x) 

{ 
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return 0.4342944819032518/x; 

} 

III Derivative of square root function, 
inline double sqrt_derivative (double x) 

{ 

return 0.5/std : :sqrt(x); 

} 

III Derivative of absolute value function, 
inline double abs_derivative (double x) 

{ 

return x > 0.0 ? 1.0 : -1.0; 

} 

} // namespace SAD 

// Add all mathematical functions to the namespace std 
// for instance, std::sin(x) also works with objects of 
namespace std 
{ 

#define IMPL_FUN_1(FUN,DER) 

inline SAD::Variable FUN(const SAD 

{ 

double val = FUN(x()); 
double der = DER(x()); 

SAD::Variable res(x); /* 

res.setValue(val); /* 

res.scaleDependencies(der); /* 

return res; 

} 


IMPL 

_FUN_ 

1 

( sin , 

SAD 

:sin_derivative) 

IMPL 

_FUN_ 

1 

(cos , 

SAD 

:cos_derivative) 

IMPL 

_FUN_ 

1 

(tan , 

SAD 

:tan_derivative) 

IMPL 

_FUN_ 

1 

(asin , 

SAD 

:asin_derivative) 

IMPL 

_FUN_ 

1 

(acos , 

SAD 

:acos_derivative) 

IMPL 

_FUN_ 

1 

(atan , 

SAD 

:atan_derivative) 

IMPL 

_FUN_ 

1 

(sinh , 

SAD 

:sinh_derivative) 

IMPL 

_FUN_ 

1 

(cosh , 

SAD 

:cosh_derivative) 

IMPL 

_FUN_ 

1 

(tanh , 

SAD 

:t anh_derivative) 

IMPL 

_FUN_ 

1 

(exp , 

SAD 

:exp_derivative) 

IMPL 

_FUN_ 

1 

(log , 

SAD 

:1og_derivative) 

IMPL 

_FUN_ 

1 

(loglO , 

SAD 

:1ogl0_derivative) 

IMPL 

_FUN_ 

1 

(sqrt , 

SAD 

:sqrt_derivative) 

IMPL 

_FUN_ 

1 

(abs , 

SAD 

:abs_derivative) 


#undef IMPL_FUN_1 
} // namespace std 


::Variablefe x) 


copy derivatives 
set function val 
compute derivati 


o that , 

type Variable. 


of X*/ 
e f(x) */ 
es of f (x) */ 


#endif // SAD_VARIABLE_0PERAT0RS_HPP 
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